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—I Abstract 

(-H In a series of recent works it has been shown that a class of simple models of evolving popu- 

I lations under selection leads to genealogical trees whose statistics are given by the Bolthausen- 

5^ Sznitman coalescent rather than by the well known Kingman coalescent in the case of neutral 

^ evolution. Here we show that when conditioning the genealogies on the speed of evolution, 

4_i one finds a one parameter family of tree statistics which interpolates between the Bolthausen- 

2 Sznitman and Kingman's coalescents. This interpolation can be calculated explicitly for one 

r1 specific version of the model, the exponential model. Numerical simulations of another version 

\_^ of the model and a phenomenological theory indicate that this one-parameter family of tree 

statistics could be universal. We compare this tree structure with those appearing in other 
contexts, in particular in the mean field theory of spin glasses. 



1 Introduction 

•^ An important question in the study of evolving populations is to understand the effect of selection 

Lj on the ancestry and on the genealogies O El HI O [1]. In absence of selection, for a well mixed 

fvq population such as in the Wright-Fisher model, the statistical properties of the genealogy of a large 

^S| population of constant size is described by Kingman's coalescent [HI [HI IS] . Recent attempts to 

f^ modify the Wright-Fisher model in order to introduce selection lead to a change of the statistical 

C^ properties of genealogies: in [TUl fTTj . the study of a whole class of models indicates that the 

T-H genealogies of populations evolving under selection are given by Bolthausen-Sznitman's coalescent 

^"^ [12] rather than by Kingman's (this has been shown analytically only for one specific version of the 

^ model, the exponential model (see below), but it has also been checked in numerical simulations 
and proved for a modified version of the model where the effect of selection is represented by a 

^ moving absorbing wall along the fitness axis [H]). In the present paper, we further study these 

5^ simple models of evolution with selection and we calculate how the statistical properties of the 
genealogies are correlated to the speed of evolution. 
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The models of evolution with selection we consider here have been introduced in [TOl E] (see 
also [13 [H])- They can be defined as follows: at each generation t the population consists of 
a fixed number N of individuals and each individual i is characterized by a single number Xi{t) 
representing its adaptation in the environment. So Xi{t) is the position of individual i on a fitness 
or adaptation axis (very much like in the Bak-Sneppen model [16)). This individual has several 
offspring at positions Xi{t) + ei.i(i), Xi{t) + £», 2(^)1 Xi{t) + €1^3(1), etc., where the eij{t) are random 
numbers representing the change of adaptation due to mutations between parent i and child j. The 
total number of offspring produced this way by all individuals at a given generation t exceeds N; 
the population at generation i + 1 is then obtained by keeping the N most adapted children (i.e. 
the N rightmost points along the axis) among all these offspring, see figure [I] The model is fully 
specified when the distribution of the number offspring of each individual and the distribution of 
the random shifts e are given. 
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Figure 1: Three time steps for a population of size N 
two offspring. 



4 in a model where each individual has 



After letting such a model evolve for a large number t of generations, the positions Xi{t) of the 
individuals on the adaptation axis from a cloud of points grouped around a position Xt which grows 
linearly with time with some velocity v^. There is some arbitrariness in the way this position Xt 
can be defined (one could choose for example Xt to be the position of the rightmost individual, or 
of the leftmost individual or the center of mass of the population) but, as the positions of all the 
individuals remain grouped, a change of definition modifies the value of Xt by an amount which 
does not grow with time and therefore does not affect the velocity vn- In addition to the velocity, 
the position Xt has fluctuations: in particular it diffuses with a variance which grows linearly in 
time [I7|. 

At each generation, one can study the genealogy of the population of this model by considering 
the matrix Tij{t) of the ages of the most recent common ancestors, or coalescence times, of all pairs 
of individuals i and j living at generation t; see figure [2] As the genealogy is a tree, the whole 
ancestry of the population can be deduced from the knowledge of this matrix. In particular the age 
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Figure 2: Left: an exemple of the genealogical tree of A^ = 5 individuals. Right: the coalescence 
times Tij{t) corresponding to the tree. 



of the most recent common ancestor of any subset of the population can be expressed in terms of 
this matrix: for example the age Ti^^kit) of the most recent common ancestor of three individuals 
i, j and k at generation t is simply 



n,3,k{t) = ra-AX {Tij{t),Ti^k{t),Tj^k{t)) ■ 



(1) 



One property common to the models of evolution under selection studied here and to the neutral 
models of evolution described by Kingman's coalescent is that the heights and the shapes of the 
genealogical trees fluctuate with t even when the size N of the population becomes large. The 
statistical properties of these trees and their time scales are however different: for instance, the 
typical age of the most recent common ancestors of k individuals grows logarithmically with the 
population size N in presence of selection (as in the models studied here) while it grows linearly 
with N in the neutral case. 

There are several ways of describing the statistical properties of these trees (see section l5| . In 
[10l|TT| we chose to characterize them by the average coalescence times (T^) of k individuals chosen 
at random in the population: 

{n)^{n,,...,At))- (2) 

(here (•) denotes an average over the individuals zi, . . . , zj. and over the generation t). The TV depen- 
dence of (T2) gives the time scale over which coalescence events occur, while the ratios {Tk) / {T2) 
are a signature of the statistical properties of the shape of the trees. We found in (TUl E] that 
these ratios in presence of selection converge when A^ — > 00 to those of a Bolthausen-Snitzmann 
coalescent: 



{T2 



5 
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(Bolthausen-Snitzmann) 



in contrast to the neutral case where they converge to those of Kingman's coalescent: 

(T3) _ 4 (T4) _ 3 



{T2 



{T2 



(Kingman) 



(3) 



(4) 



Our goal here is to calculate how these ratios are correlated to the speed of evolution, by 
weighting all the events during a long time interval r by a factor e~°'^^ . {^ < favors events with 



a speed of evolution faster than average, while /3 > correponds to events with a slower speed 
of evolution.) Our main result, derived below for the exponential model, is that the above ratios 
become for large N 



in) 5 + 4/3 (T4) 100 + 204;3 + 133;32 ^_ 27/33 



(Ta) 4 + 3;9' (T2) 72 + 142/3 + 90^2.^18/33 



(5) 



It is remarkable that these expressions interpolate between the neutral case (Kingman) for /3 — >■ +00 
(low speed limit) and the selection case (Bolthausen-Snitzmann) for /3 = 0. When /3 ^- — 1 (high 
speed limit), all the ratios become 1 indicating a "star-shaped" coalescent. 

This paper is organized as follows: in section [2j we explain how the weighting by the factor 
g-^^T is done. In section pi we show that in presence of the bias, one version of the model (the 
exponential model) can be solved exactly by analyzing a coalescent model, the rates of which depend 
on /3. This leads to ([5|. In section |4] we argue using the phenomenological theory developed in [11] 
that (l5| should remain valid for other versions of the model up to a change of scale of /3. Lastly 
in section [5] we compare the /3-random tree structure which leads to ((sl) to the statistics of the 
partitions in mean field spin glasses and in the Poisson-Dirichlet distribution. 

2 How to condition on the velocity 

If one performs a simulation of the model described in the introduction, one can measure at each 
generation t the position Xt of the population (defined in any reasonable way: as explained in 
the introduction, the precise definition does not matter) and the ages T2(i), . . . ,Tk{t) of the most 
recent common ancestor of 2, . . . , fc individuals chosen at random in the population at time t (for 
more efficiency one can average these times Tk{t) over all the choices of the k individuals in the 
population at time t). Then we choose a long time interval r and we want to determine 

(^^)^ = J™.-L ^,-px.^ ■ (6) 



t=l 



2.1 Theoretical considerations 

As the model is a Markov process and correlations decay fast enough in time, we expect that, for 
large t, 



The knowledge of G(/3) determines all the cumulants of the position Xt 

(8) 



Ih^ i^^ ^ (_)« d"G(/3) 



/3=0 



t ' ' d/3« 

It is also related to the large deviation function F{v) of the velocity defined by 

Proba(Xt = vt) - e*-^'") (9) 

through a Legendre transform 

G(/3) = max[-^w + F{v)]. (10) 



For large t, the value of v which dominates the weighted averages in ^ and ([7| is given by 

with fluctuations of order i^^/^. Therefore the weighted averages m» become equivalent in the 



i — >■ cxD limit to conditioning on the velocity v given by (11). This is, in the present context, the 



analog of the well known equivalence of ensembles in statistical physics. 

2.2 In numerical simulations 

Numerically it is difficult to perform averages such as ^ because the events which dominate both 
the numerator and the denominator of (pi) are rare events. In order to overcome this difficulty, we use 
an importance sampling method. We consider a sample periodic in time of period t where r is chosen 
large enough. This means that the random shifts eij(i) are periodic in time {£i,j{t + T) — ei,j{t) for 
all i, j and t). With these periodic eij{t), the evolution of the system becomes also periodic in time: 
the shift Xt of the position of the population after one period r can therefore be unambigously 
defined and depends on all the eij(t). Then, we perform a standard Monte-Carlo simulation: at 
each step we try a new sample by changing some of the eij(t) and we let the system evolve till it 
becomes periodic (here we change all the eij{t) at a random time t uniformly distributed between 
1 and r). The outcome of this change is to modify the Xr to a new value X"°^. Then, as always 
with a Metropohs algorithm, we accept the change with a probability max [l, exp[— /3(X"''^ — X^.)]] . 
With this procedure samples are produced with a weight exp[— /3Xt], so that by averaging quantities 
such as the Tk over many samples one gets an estimate of ( 6 ) . 

We have simulated the exponential model (see section 3 where we give the precise definition 
of the exponential model and its analytic solution in the iV — >■ oo limit) for N = 100 and a value 
of r sa 301nA^ which is much larger than {T2) [11]. For each value of /3, we measured (T2), {T3), 
(T4) averaged on 10^ Monte-Carlo steps. We have also simulated a more generic model where each 
individual has exactly two offspring with independent random shifts ctjit) uniformly distributed 
between and 1. This model cannot be solved exactly, but a phenomenological theory (see [TT] 
and section 111) predicts when N ^ 00 the same statistics of the genealogical trees ^ as in the 
exponential model with /3 replaced by /3/7, where 7 ~ 5.262 is the value which minimizes the 
function ln[2(e^ - l)/7]/7, see section |4J We simulated the sizes iV = 30, iV = 100 and N = 300 
with values of r ~ 8 In^ N which is much larger than (T2), and again averaged over 10^ Monte-Carlo 
steps. We also checked for several values of /3 that our results remain unchanged by choosing a value 
of the period r twice as big (results not shown), indicating that our Monte-Carlo results would be 
the same for an infinite time-window. 

The results for {T^) / {T2) and {T^) / {T2) are presented in figures Is] and El As in [TT], we observe 
that for the exponential model, the results are already very close to their asymptotic limits even 
for N = 100; only for /? < — 1 there is a discrepancy with the theoretical prediction ([5]). As in [TT) . 
the convergence is however slower in the generic case but the curves seem in both cases to converge 
to the prediction. 

3 Exponential model 

In this section we consider a version of the model, the exponential model, which can be solved 
exactly [TT]. In the exponential model, the shifts of the offspring of each individual are generated 
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Figure 3: [T^) /{T2) as a function of /? for the exponential model and the prediction (l5|, and as a 
function of /3/7 for the generic model described in the text. 



by a Poisson process of density p(e) — e~'^ . This means that an individual at position Xi{t) has a 
probability e~^de of having an offspring in the infinitesimal interval {xi{t) + e, Xi{t) + e + de). Then 
the population at the next generation is obtained by selecting the N rightmost points among all the 
offspring produced by generation t. (Note that in the exponential model the number of offspring 
produced by each generation is infinite but their number at the right of any position y is finite. 
There is therefore no problem to select the A'' survivors at generation t + 1). 

In the exponential model, there is a convenient way of defining the position Xt of the population 



X, =ln 



N 



,'^^it) 



(12) 



The simplicity of the exponential model comes from the fact that, with this definition of Xt, one 
has 



E' 



-[x-xdt)\ 



,-(^-^t) 



(13) 



which means that one can generate the offspring of the whole population at time t by replacing the 
N Poisson processes centered at the positions Xi (i) by a single Poisson process centered at position 
Xt. Therefore, with definition (12 1 of Xt, the N points Xi{t + 1) at generation i + 1 are the N 
rightmost points of a Poisson point process with a density exp[— (x — Xt)]. 

As explained in [111 a way of drawing these N points is to choose a number z with a density 
of probability Proba(2:) = exp[— (A^ + l)z — e~^]/N\ and, independently, N numbers yi with an 
exponential density Proba(y) = e^y0{y); the points Xi{t + 1) are then given (in an arbitrary order) 

by 

x,{t + l)=Xt + z + y„ (14) 
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Figure 4: {T^} /{T2) as a function of /? for the exponential model and the prediction (l5|, and as a 
function of /3/7 for the generic model described in the text. 



and one gets from (12) 



r N 



Xt+i = Xf + z + In 



Y^ey 



(15) 



We see that, with the definition (12 1 of Xt, the differences Xt+i — Xt are independent variables. 
Therefore from 



and G{P) can be computed by averaging over a single generation 



eG(/3) = (e-^- 



AT 



^e- 



1=1 



r(iv + i + /3) 



r(7V + i 

Using the integral representation (valid for /3 > 0) 

1 



dyi ■ ■ ■ / dj/„ e 
JQ 



-Vl Vn 



N 



^e- 



r(/3) Jo 



A-P^^^ dAA^-^e-^-^ 



(for /3 > 0), 



one obtains 



r(iv + i)r(/3) io 



w 



(17) 



(18) 
(19) 



where the integral Io{X) and more general integrals Ip{X) are defined as 



Ip{X)= / dye'^P-^^y^^"" ^X^^P 

Jo J\ 



duuP-^e-"". 



(20) 



(These integrals are in fact, up to a simple change of variables, incomplete gamma functions). 
For large N, the expression (19) is dominated by small values of A where the integrals Ip{X) for 
non-negative integers p can be approximated by |llj : 



/o(A) = 1 + A(ln A + 7B - 1) + O(A'), /i(A) = -(In A + 7^) + 0(A), 
Ip>2{X)=^^^+0{X'-P), 



(21) 



xp- 



where je — ^r'(l) ~ .577 is Euler's constant. 



Given the small A expansion of /o(A), the integral in (191 is dominated, for N large and /3 of 

/{N In N] 

/o(A)^ - e-^ 



order 1, by A of order l/(iVln A^). Making the change of variable A = fi/{N\nN) one has 

In/i — Inln A^ + 7£; — 1 ///ln/;\2 



1 + Ai- 



InA^ 



VlnA^ / 



22 



We can now evaluate ^ for N large and /3 of order 1; using ([22]) and TjN + l + /3)/r{N + 1) ~ 7V^, 
one gets [TT] 



=G(/3) _ 



1 



1 



/3 /r'(/? + i) 



or 



G(/3) 



In'' TV 
-/31nlniV- 



InA^ \^r(/3 + l) 

f3 /r'(/3 + i) 



In TV Vr(/3 + i) 



- In In A^ + 7_E - 1 
In In A^ + 7_E - 1 



(23) 



(24) 



-G'{(3), see (11), varying /3 does not change the leading A^ dependence 



We see that, as v 

w ~ In In A^ of the velocity but only shifts it by a small amount of order In In N/ In A^ which vanishes 
in the N -^ 00 limit. We are now going to show that, on the contrary, (3 does change the statistical 
properties of the trees even in the A^ — > 00 limit. 



3.1 Trees 

We have already seen that all the offspring produced by generation t are distributed according to 
a Poisson point process of density exp[— (a: — Xt)]. On the other hand, the offspring of individual 
Xi{t) are distributed as a Poisson point process of density exp \—[x — Xi{t)]\. This implies that, 
given that there is an offspring in an interval dx around x, its probability of being an offspring of 
Xi{t) is 



W, 



,x,it)-Xt ^ 



^x,{t) 



^- ,.,(*) 



(25) 



This probability is independent of x. Therefore the probability Qp{t) that p individuals at generation 
t + 1 have the same ancestor at generation t is 



AT 



JV 



gp(t) = ^W^f = ^eP-'W-P^'. 



(26) 



If one weights these coalescence rates with the factor e ^^'^ , then using ( 14 ) and (151 with t replaced 



by i — 1 and using the fact that Xt-i, z and the yi are independent, one gets 

/,„(n-(«+rtXA /^m-(»+rlln[E.4.e..]\ 



with the yi independent exponential variables. The numerator and the denominator can be com- 



puted in the same way as in ( 19 1 and one obtains 



^Q^^^-^rWTp) J,-'dxx^-Ho{xr ■ ^ ^ 

We take p > 2 and /? of order 1. The integrals are dominated by A = ^/{NlnN) and /j, of order 1. 

To the leading order, Io{X)^ ~ e~^, see ( [22| , and using (21) for Ip{X) one easily gets, to leading 

order 

^ I (p-2)!r(/3 + l) ^ 1 (p-2)! 

^^P^^-lnAT r(/3 + p) IniV (l + /3)(2 + /3)---(p-l + /3)- ^ ' 

After rescaling time by a factor InA^, one gets a coalescent with transition rates qp — {Qp)i3\nN. 
One can check that 

qp = ^^~p?!^^^N^^^ = / xP-^Hdx) with A(da;) = (1 ^ x)^ dx, (30) 



Using the expressions ( [55| ) of the appendix, where the ratios {T^i) / {T2) and {T4} / {T2) have been 
obtained for a general coalescent, one finally gets ([5|. 

4 The phenomenological theory 

In this section, we show that the phenomenological theory developed in [17l [11] in the context of the 
noisy Fisher-KPP equation predicts that ([5| remains valid for other versions of the model described 
in the introduction. When the number of offspring of each individual is bounded and when the shifts 
eij are also bounded, one can describe the evolution of the population by a noisy traveling wave 
equation of the Fisher-KPP type. In [IZIITT], a phenomenological theory was proposed to describe 
the large N behavior of these noisy equations. In the N —>■ 00 limit, the effect of noise vanishes and 
the traveling wave has a finite velocity Vao (in contrast to the exponential model where the velocity 
diverges as A^ — >■ 00). The first correction when N is large can be understood by considering the 
cutoff introduced by the discrete number of particles: this leads to Wcutoff = Voo — A/ In N. The 
next order correction leads to a positive term of order In In A/ In A which can be understood, as 
well as the fluctuations of Xt, by the following phenomenological theory: the front has, most of 
the time, the shape and the velocity predicted by the cutoff theory. However, every (typically) 
In'^ A time steps, a rare event occurs where some particles escape significantly ahead of the front. 
When this happens, the shape of the front is at first deformed, but it relaxes to its cutoff shape 
after ~ In A time steps. The end result is a finite increase of the position of the front [TT]. It has 
been shown that this phenomenological theory predicts genealogies described by the Bolthausen- 
Snitzmann coalescent. We are now going to show that when we condition on the velocity by using 
the weight exp(— /^X^), this leads to ([s]). 



We consider a time interval At which is large compared to In N but small compared to In N. 
During this interval, there is a small probability p{S) dS At that an event of size 6 occurs. When 
this happens, the front position increases (after relaxation) by R{S). The time interval In N <^ 
At <C In'^ N is such that each event has the time to relax during At and that the probability that 
two events occur during the same time interval is negligible. 

With these notations, the position Xt of the front evolves according to 



X, 



t+At ■ 



X, 



^cutoff Af + R{S) proba. p{S) dS At, 
fcutoff Ai + proba. 1 — Atjd5p(5). 



We argued in [TT] that, for large 6, 

p{6) « Cle-''^ R{6) 



In 1 + Q 






CiC2« Af"(7). 



(31) 



(32) 



The number 7 and the function v{r) depend on the details of the model. (For N —^ 00, v{r) gives 
the velocity of a front starting with the initial condition e^''^. For a step initial condition, the 
system moves at the velocity ^(7) where 7 is the value at which v(r) reaches its minimum.) 
Using (31| to compute G(/3) given by ([7|, one gets 



G(/3) = -/3z;cutoff + / d5p{5) [e-^«(^) - 1 



Let us now weight all the events by the factor exp(~/3Xj). (31 1 becomes 



X.. 



t+At 



Xt^ 



WcutoffAi 
Wcutoff Ai 



R{5) proba. 




proba. ^e-''^-'""'^* [l " At / dSp{S)] . 



(33) 



(34) 



Z{f3) is such that the probabilities are normalized; clearly Z{/3) = e^*'^^^\ 

Now, we can try to determine the probability Qp At that the p particles coalesce into one 
during the time interval At. We argued in [111 that when a rare event of size 6 occurs, a fraction 
/ = l — e~^^-^^'> of the population is replaced by the offspring of the single particle that originated the 



event (7 is the model specific number appearing in (32)). When this happens, there is a probability 



/'' that the p particles coalesce during that interval of time At. This leads to 



tt/p 



dSp{6)e-^'"^'^ \l - e-^«(*) 



We can now use (32) in (35). Rewriting the integral in term of the variable / = 1 — e 



(35) 



-7i?,(5) 



(the fraction of the population replaced by the offspring of an individual) , one gets p{S) d6 
CiC2-f^/ln^Ndf/p, so that 



In^iV 



F-'{i~ff/-'df, 



(36) 



which is the same as ( 30 ) up to a prefactor (which only changes the time scale) and the fact that 



/3 is replace by f3/"f. Therefore the phenomenological theory leads to the same coalescent as in the 
exponential model but with a different time scale of order In'^ N instead of In TV. 



10 



5 Comparison with mean-field spin glasses and the Poisson- 
Dirichlet distribution 

5.1 Various ways of characterizing random trees 

There are several ways of characterizing the statistical properties of the trees generated by some 
given coalescence rates. (We only consider here the cases where the coalescence rates do not vary 
in time, where the particles play symmetrical roles and where at most one coalescence event can 
occur during an interval of time dt.) 



One can specify the coalescence rates qp (which take the values ( 30 1 for the models of evolution 
with selection that we consider in this paper). In terms of these coalescence rates, Kingman's 
coalescent corresponds to 

92^0, qp^O forp>3, (37) 

while the Bolthausen-Sznitman coalescent corresponds to 

qp = ^. (38) 



It is easy to see that the rates ( 30 ) interpolate between ( 37 ) for /3 = cx) and ( 38 ) for P — 



• One can alternatively specify all the ratios {Tp) /{T2) ■ It is clear (see ( 52|55 ) in the appendix) 
that the knowledge of the qp determines all these time ratios and conversely that the knowledge 
of the time ratios allows one to calculate all the ratios qp/q2- 

• One can also characterize the trees by the partition of the population they induce at a given 
time in the past: the population at a generation t can be decomposed into several r-families 
where, by definition of these families, two individuals i and j belong to the same r-family if 
the age of their most recent common ancestor is less than r (i.e. Tij{t) < t). One can then 
associate to this r-partition of the population at generation t the following numbers 

Yt\t)^{e{r-r,,_,,))^ (39) 

where {■)t means an average over all the possible choices of the k individuals ii,...,ik at 
generation t. 

One can interpret these Y^'^ (t) as the probability that k individuals chosen at random in 
the population at generation t belong to the same r-family. These i'fc (i) fluctuate from 

generation to generation and the expressions {Y^'^ ) of their averages over t can be computed 
in terms of the coalescence rates qp. In the appendix, they are given for /c = 2, 3,4 by the 
quantities Zfc^i(T) = (Y^^'). 

Knowing all the (Y^"^ ) (even for a single t) determines also in principle all the coalescence 
rates and therefore all the statistical properties of the trees. 

5.2 Comparison with mean-field spin glasses 

Very much like in the coalescence problems discussed above, where all the individuals at a generation 
t can be grouped into r-families, one can group the spin configurations of a spin glass model 

11 



according to their distances d (or to their overlap= 1 — d) in phase space. One can then define 
[TS], for a given sample, the probabihty Yk(d) that k configurations, at thermal equilibrium, have 
all their k{k — l)/2 mutual distances in phase space less than d. 

One of the predictions [TSlIinillS] of the Parisi solution [53J|2TJ|21] of the Sherrington Kirkpatrick 
model [551 Hi] is that this Yfe(d) fluctuates with the spin glass sample even when the system size 
becomes large. The Parisi theory predicts also all the statistical properties of these Yk{d). For 
example 

(n(.))^lim, ,(.,,). ^.^^^^ (40) 

where according to the broken replica symmetry 

r(l - n) r(fc - fi) 

^^^"'^^ = r(fc-n)r(i-M) - ^''^ 



In (40) all the dependence on the distance d, on the details of the model, and on the parameters 
such as the temperature or the magnetic field is through the parameter /i. Formula ( [4l] ) follows 
from a very simple replica calculation: assume that one has n replicas grouped into n//i families of 
jU replicas, yk{n, //) is simply the probability that k replicas chosen at random among the n replicas 
belong to the same family. 



For typical samples one has to take the n — >■ limit as in ( 40 1 and the statistics of the Yfe coincide 
with those of the Bolthausen-Sznitman coalescent: one can check that the Zk^i{T) = (^^^ ) 



obtained in (57) coincides with (40) by choosing /i = e"''^'^ and the qp given by (38). 

In the spin glass case, one can also weight the samples according to their free energy (by 
weighting them by a factor Z" where Z is the partition function [27l[26]). One then expects from 
the replica theory [26 that the statistics of the Yfe (d) be modified and that 

^-'^m^-V.M. (42) 



One can check easily from ( 57 ) with qp given by ( 30 ) that there exists no choice of n and /i as 
functions of /3 and t such that yk{n,^) = Zk^iir). Therefore, although the statistical properties 
of the trees in the model of evolution with selection and in the spin glass problem are the same for 
typical samples, they become different when one introduces a bias (related to the free energy in the 
spin glass problem and to the speed of adaptation in the models of evolution with selection) . 

5.3 The Poisson-Dirichlet distribution 

The Poisson-Dirichlet distribution [5SJ [55] is a probability distribution of the partitions of a unit 
interval into infinitely many subintervals. It is parametrized by two parameters a and 0. One way 
of defining the Poisson-Dirichlet distribution is to consider an infinite sequence zi, 2:2, . . . , Znj ■ • ■ of 
independent numbers, each z„ being distributed according to a distribution 

I (1 — a)i [0 + na) 

(which is a /3 distribution) . Then one considers a partition of the unit interval into subintervals of 
lengths Wi, W2, . . . Wn, ■ ■ ■ with 

Wi + W2 + --- + Wn + --- = l, (44) 
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where the Wi are given by 

(45) 

Wn = (1- Zi){l - Z2) • • • (1 - Zn^l)z„, 

For such a partition one can introduce the quantities 

n- = Ew^'' (46) 

i 

which represents the probabihty that k points chosen at random on the unit interval faU in the 
same subinterval. It is easy to check that when one averages over the z^, one gets 

{Yk)a.fi = (^l) + ((1 - Zl)')(n)a,a+9, (47) 

The sohition of this recursion is 

_ ril + 9)T{k-a) 



One can notice [30] that these expressions are identical to the replica expressions (41) of yfe(n,/x) 
when one chooses 9 = —n and a = fi. Therefore as soon as one introduces the bias /3 f^ the 
statistical properties r-families of our models of evolution with selection differ from those of the 
Poisson-Dirichlet distribution. 

6 Conclusion 

In this paper we have seen that, for a family of simple models of evolution under selection, the 
statistics of the genalogies are modified when conditioning on the speed of evolution ^. For one 
particular version of these models, the exponential model, the trees can be generated by a coalescent 
with modified rates ( 30 1 . Numerical simulations (figures p^ and El and a phenomenological theory 



(section HI) indicate a similar behavior of more generic versions of the model. 

Despite their simplicity, there is not yet a full theoretical understanding of the models of evolu- 
tion with selection we consider here. The introduction of the bias opens new questions which would 
be interesting to consider. For exemple, what is the effect of the bias on the steady state density 
profile of the population along the fitness axis, or on the distances between the rightmost points 
in the population? Numerically, the Monte-Carlo approach we developed here should give a rather 
powerful tool to study these questions accurately and to test more precisely the phenomenological 
theory developped in [T71 [TT] . 

It would also be interesting to study the genealogies of other models of evolution with selection [5T] 
to test the genericity of our results. 

We are very happy to dedicate this work to David Sherrington, on the occasion of his 70*^ 
birthday. 
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A Coalescence times and sizes of families in the A coalescent 

In this appendix we calculate a few simple properties of a continuous time coalescent defined as 
follows: one starts with N points and during every infinitesimal time interval di ^ 1, every subset 
of k points has a probability q^ dt of coalescing into one point. It is assumed that there is at most 
one coalescence event during a time interval dt. 

This model is called the A-coalescent [331 1321 1311 [ISl US] . The coalescence rates can be written 
[33] in terms of a positive measure A on the interval (0, 1) 



Qk 



„fc-2 



A(da;), 



(49) 



and, more generally [33], the rate Xb^k at which the k > 2 first points out of 6 coalesce into one 
point (while the other b — k points remain single) is given by 



«1 b—k 

Xb,k^ / x''~'{l-x)''-''A{dx)^Yl 



{b-k)l 



— ' n\(b — k — 72)\ 
n=o ^ ' 



(-l)"9n+fe. 



(50) 



It is more convenient in the following to use of the quantities rh(b'), defined as the rate at which a 
set of b points coalesce into a set oi b' < b points. Clearly 



6! 



""(^'^ ~ ib'-mb-b' + iy.^'-'-''-'' 



(51) 



(This simply means that the total number of distinct points jumps from b to 6' with probability 
ri,{b')dt during an infinitesimal time interval dt. The binomial factor in ( [51[ ) comes from the number 
of ways of choosing the b — b' + 1 points which coalesce.) 

As already noticed, the rates (30 1 correspond to A(da;) = (1 — x)^dx. 



By analyzing what happens during a time interval dt one can then see that the age Ti, of the 
most recent common ancestor of b individuals chosen at generation t + dt satisfies 



Tbit + dt)= < 



dt 

dt + r2(t) 

dt + rb_i(t) 
I dt + rb(t) 



with probability rb(l)dt, 
rb(2)dt. 



rb{b-l)dt, 



(52) 



Therefore one can determine recursively the average coalescence times by writing that {Tb{t + dt)) = 
{Tb{t)), which leads to 



b-l 



E ^»(b') 



.b' = l 



6-1 



{Tb) = l + J2rbib'){Tb,) 



6' =2 



As (50 511 imply that 



''2(1) = 92, 
''3(1) = 93, 
''4(1) = 94, 



r3(2)=3((Z2-93), 
r4(2)=4(g3-(j4), 



7-4(3) = 6(92 -293 +94), 



(53) 



(54) 
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one gets 

^y^^ ^1 Os) ^ 4g2 - 3g3 

92' (72) 892-253' 

(r4) ^ 27g| - 56g2g3 + 28gi + 12g2g4 - 10g3g4 
(T2) (3g2-2g3)(6g2-8g3 + 3q4) 

If one defines Zb^hi (t) as the probabifity that b points have coalesced into b' points during some 
time T, one can easily see that it evolves according to 

dZb^b' ^ 
dr 

b">b' b"<b 



J2 n"{h') Z,^y, ~ Y, n'{b") Zk^,,, (56) 



with the initial condition that Zt,_j.ft'(0) = (Jb.t,'. Then, using the expressions (54) one gets 



^2^2 

Z3^2 = -e-"^" - -e-(3?2-2,3)r 

^3^1 = 1 - ^e-'^^^ + le-(3«^-2«^)", (57) 

7 _ „-(6g2-8g3+3q4)T 

•^4-5-4 — t; , 

^^^^ ^ 9g2-14g3 + 5g4 ^-,,. _ 2g-(3g.-2,3)r ^ 6g2 - 10g3 + 4g4 ^,(6q,-8q3+3<74)r ^ 

5^2 - 8^3 + 3^4 5^2 - 8^3 + 3^4 

Z^^^ = 1 _ 9g2-14g3+5g4 ^-g,^ _^ ^_(3q,_293)r _ g2 " 2^3 + 94 ^-(6^,-8a3+3q4)r ^ 

^ 5^2 - 8(73 + 3(?4 5^2 - 8q'3 + 3g4 
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